/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/*!
    \file umfmatrix.hpp
    \brief This file contains the UMFMatrix class
 *
 *  This file contains the UMFMatrix class
 */


#ifndef UMFMATRIX_HPP
#define UMFMATRIX_HPP 1

#include<valarray>
#include<fstream>

extern "C" {
    #include<umfpack.h>
}

namespace spectral {

/**
 * \class UMFMatrix
 * \brief Store the UMF matrix
 *
 * CSR (Compact Sparse Row) Matrix format and UMFPACK wrappers
 */

/**
 * \fn UMFMatrix::UMFMatrix( )
 * \brief Default constructor
 */

/**
 * \fn UMFMatrix::UMFMatrix(int sizeA,
 *                          int sizeNZ)
 * \brief The UMFMatrix class constructor
   \param sizeA Matrix dimension
   \param sizeNZ Number of non zero elements
*/

/**
 * \fn UMFMatrix::UMFMatrix(UMFMatrix const &)
 * \brief The UMFMatrix class copy constructor
   \param Class to be copied
*/

/**
 * \fn UMFMatrix::UMFMatrix & operator=(UMFMatrix & B)
 * \brief The overloading of the = operator for UMFMatrix class
   \param B Class to be copied
*/

/**
 * \fn ~UMFMatrix()
 * \brief Destructor
 */

/**
 * \fn void UMFMatrix::resize_sizeA(int sizeAnew)
 * \brief Resizes the dimension of the problem
 * \param sizeAnew The size of the new problem
 */

/**
 * \fn void UMFMatrix::resize_sizeNZ(int sizeNZnew)
 * \brief Resizes the number of non zero elements of the problem
 * \param sizeNZnew The number of non zero elements of the problem
 */

/**
 * \fn void UMFMatrix::resize(int sizeA,int sizeNZ)
 * \brief Resizes the problem
 * \param sizeAnew The size of the new problem
 * \param sizeNZnew The number of non zero elements of the problem
 */

/**
 * \fn void UMFMatrix::solve(std::valarray<double> const & b,
                             std::valarray<double> & u) const
 * \brief Solves the problem
 * \param b The known term vector
 * \param u The solution vector
 */

/**
 * \fn void UMFMatrix::solve(double const * b,double * x) const
 * \brief Solves the problem
 * \param b The known term vector
 * \param x The solution vector
 */

/**
 * \fn int UMFMatrix::getsizeA() const
 * \brief Returns sizeA
 */

/**
 * \fn int UMFMatrix::getsizeNZ() const
 * \brief Returns sizeNZ
 */

/**
 * \fn int &  UMFMatrix::getAp(int idxpos)
 * \brief Gets the element of Ap
 * \param idxpos The index of the element
 */

/**
 * \fn int &  UMFMatrix::getAi(int idxpos)
 * \brief Gets the element of Ai
 * \param idxpos The index of the element
 */

/**
 * \fn int &  UMFMatrix::getAx(int idxpos)
 * \brief Gets the element of Ax
 * \param idxpos The index of the element
 */

/**
 * \fn void UMFMatrix::multiply(std::valarray<double> const & in,
                                std::valarray<double> & out) const
 * \brief Multiplies a vector for the matrix
 * \param in The vector that will be multiplied
 * \param out The solution
 */

/**
 * \fn void UMFMatrix::multiply(double const * V, double * AV) const
 * \brief Multiplies a vector for the matrix
 * \param V The vector that will be multiplied
 * \param AV The solution
 */

/**
 * \fn UMFMatrix::UMFMatrix & operator+=(UMFMatrix & B)
 * \brief The overloading of the += operator for UMFMatrix class
   \param B Matrix to be added
*/

/**
 * \fn UMFMatrix::UMFMatrix & operator*=(double const & alpha)
 * \brief The overloading of the *= operator for UMFMatrix class
   \param alpha Constant to be multiplied
*/

/**
 * \fn UMFMatrix::UMFMatrix & operator+(UMFMatrix A, UMFMatrix B)
 * \brief The overloading of the + operator between two UMFMatrix objects.
 */

/**
 * \fn UMFMatrix::UMFMatrix & operator*(UMFMatrix A, double alpha)
 * \brief The overloading of the * operator for the multiplication
 *        of a UMFMatrix with a double.
 */

/**
 * \fn UMFMatrix::UMFMatrix & operator*(double alpha, UMFMatrix A)
 * \brief The overloading of the * operator for the multiplication
 *        of a UMFMatrix with a double.
 */

/**
 * \fn void UMFMatrix::printPattern(char const * name) const
 * \brief Prints the pattern of the matrix as a string of char
 * \param name The name of the file storing tha pattern
 */

/**
 * \fn double * * UMFMatrix::fullmatrix() const
 * \brief Copies the UMF matrix in an array of arrays: use only for debugging purposes!
 */

/**
 * \fn void UMFMatrix::factorize()
 * \brief Factorizes the matrix
 */

/**
 * \fn void UMFMatrix::clean()
 * \brief Cleans UMFPACK data
 */

// This stores the transposed matrix!!!
class UMFMatrix {
    public:
        UMFMatrix();
        UMFMatrix(int sizeA, int sizeNZ);
        UMFMatrix(UMFMatrix const &);
	UMFMatrix & operator=(UMFMatrix & B);

        ~UMFMatrix();

	void resize_sizeA(int sizeAnew);
	void resize_sizeNZ(int sizeNZnew);
	void resize(int sizeA,int sizeNZ);

	void solve(std::valarray<double> const & b,
                   std::valarray<double> & u) const;

        void solve(double const * b,double * x) const;

        int getsizeA() const;
        int getsizeNZ() const;

        inline int & getAp(int idxpos);
	inline int & getAi(int idxpos);
	inline double & getAx(int idxpos);

	void multiply(std::valarray<double> const & in,
		      std::valarray<double> & out) const;
	void multiply(double const * V, double * AV) const;

	UMFMatrix & operator+=(UMFMatrix & B);
	UMFMatrix & operator*=(double const & alpha);

	void printPattern(char const * name) const;

        // Use only for debugging purposes!
        double * * fullmatrix() const;

        void factorize();

        // Clean UMFPACK data (to be used when changing the Matrix!)
        void clean();
    private:
        // Matrix data
        int* _Ap;             /**< \brief Array of the indexes of the first
                                          element of each row */
        int* _Ai;             /**< \brief Array of the indexes of the horizontal
                                          position of the element */
        double* _Ax;          /**< \brief Array of the value in the subproblem
                                          matrix */

        // Matrix size
        int _sizeA;           /**< \brief Dimension of the problem */
        int _sizeNZ;          /**< \brief Number of non zero elements */

        // UMFPACK data
        void * _Symbolic;     /**< \brief UMF internal variable */
        void * _Numeric;      /**< \brief UMF internal variable */

        mutable double * _Control;     /**< \brief UMF control parameter */
        mutable double * _Info;        /**< \brief UMF control parameter */

        mutable int _status;  /**< \brief UMF internal variable */
        mutable int _sys;     /**< \brief UMF internal variable */

        bool _isFactorized;   /**< \brief Flag for the matrix factorization */
};

// ====================================================
// End of proper header
// ====================================================

inline int & UMFMatrix::getAp( int idxpos )  {
    assert(idxpos>=0 && idxpos<=_sizeA);
    return _Ap[idxpos];
}

inline int & UMFMatrix::getAi( int idxpos )  {
    assert(idxpos>=0 && idxpos<_sizeNZ);
    return _Ai[idxpos];
}

inline double & UMFMatrix::getAx( int idxpos )  {
    assert(idxpos>=0 && idxpos<_sizeNZ);
    return _Ax[idxpos];
}

UMFMatrix & operator+(UMFMatrix A, UMFMatrix B);

UMFMatrix & operator*(UMFMatrix A, double alpha);
UMFMatrix & operator*(double alpha, UMFMatrix A);

} // namespace

#endif
